knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(cache = TRUE)
library(tidyverse)
library(cowplot)
library(broom)
library(readxl)
library(janitor)
library(plotrix)
library(here)
library(growthTools)
library(rootSolve)
Read in data
treatments <- read_excel(here("data-general", "ChlamEE_Treatments_JB.xlsx")) %>%
clean_names() %>%
mutate(treatment = ifelse(is.na(treatment), "none", treatment)) %>%
filter(population != "cc1629")
nitrate <- read_csv(here("data-processed", "nitrate-abundances-processed.csv"))
this is the step that gets us the growth rate estimates
growth_rates_n_AICc <- nitrate %>%
filter(population != "COMBO") %>%
mutate(ln.fluor = log(RFU)) %>%
group_by(well_plate) %>%
do(grs = get.growth.rate(x = .$days, y = .$ln.fluor,id = .$well_plate, plot.best.Q = F))
Pull out the things we want
growth_sum_n_AICc <- growth_rates_n_AICc %>%
summarise(well_plate, mu = grs$best.slope, n_obs = grs$best.model.slope.n,
slope_r2 = grs$best.model.slope.r2,
best.model_r2 = grs$best.model.rsqr,
best.model = grs$best.model, best.se = grs$best.se,
contents = grs$best.model.contents)
key <- nitrate %>%
select(well_plate, population, nitrate_concentration) %>%
distinct(well_plate, .keep_all = TRUE)
AICc_growth_rates <- growth_sum_n_AICc %>%
select(-contents) %>%
mutate(IC_method = "AICc")
AICc_growth_rates2 <- left_join(AICc_growth_rates, key, by = "well_plate")
# write_csv(AICc_growth_rates2, here("data-processed", "nitrate_exp_growth_w_growthtools_AICc.csv"))
exp_params <- growth_sum_n_AICc %>%
unnest(contents %>% map(tidy, .id = "number")) ## pull out the slopes and intercepts etc.
exp_params_aug <- growth_sum_n_AICc %>%
unnest(contents %>% map(augment, .id = "number")) %>% ### now add the fitted values
rename(days = x)
exp_wide <- exp_params %>%
spread(key = term, value = estimate)
all_preds <- left_join(exp_params_aug, exp_wide)
all_preds2 <- left_join(all_preds, key, by = "well_plate") %>%
group_by(well_plate) %>%
mutate(B1 = mean(B1, na.rm = TRUE)) %>%
mutate(B2 = mean(B2, na.rm = TRUE)) %>% ### here we do some wrangling to get the colour coding right for our plots
mutate(exponential = case_when(best.model == "gr" ~ "yes",
best.model == "gr.sat" & days < B1 ~ "yes",
best.model == "gr.lag" & days < B1 ~ "yes",
best.model == "gr.lagsat" & days < B2 & days > B1 ~ "yes",
TRUE ~ "no"))
all_preds2 %>%
ggplot(aes(x = days, y = .fitted, group = well_plate, color = best.model)) + geom_line() +
geom_point(aes(x = days, y = y, shape = exponential)) +
facet_grid(nitrate_concentration ~ population) + ylab("Ln(RFU)") +xlab("Days")
Fit the Monod model to the growth rates
monod_fits <- AICc_growth_rates2 %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
rename(estimate = mu) %>%
group_by(population) %>%
do(tidy(nls(estimate ~ umax* (nitrate_concentration / (ks+ nitrate_concentration)),
data= ., start=list(ks = 1, umax = 1), algorithm="port", lower=list(c=0.01, d=0),
control = nls.control(maxiter=500, minFactor=1/204800000))))
get the fitted values
prediction_function <- function(df) {
monodcurve<-function(x){
growth_rate<- (df$umax[[1]] * (x / (df$ks[[1]] +x)))
growth_rate}
pred <- function(x) {
y <- monodcurve(x)
}
x <- seq(0, 1000, by = 0.1)
preds <- sapply(x, pred)
preds <- data.frame(x, preds) %>%
rename(nitrate_concentration.x = x,
growth_rate = preds)
}
bs_split <- monod_fits %>%
select(population, term, estimate) %>%
dplyr::ungroup() %>%
spread(key = term, value = estimate) %>%
split(.$population)
all_preds_n <- bs_split %>% ### here we just use the fitted parameters from the Monod to get the predicted values
map_df(prediction_function, .id = "population")
all_predsn_2 <- left_join(all_preds_n, treatments, by = c("population")) %>%
distinct(ancestor_id, treatment, nitrate_concentration.x, .keep_all = TRUE)
all_growth_n2 <- left_join(AICc_growth_rates2, treatments)
all_growth_n2 %>%
mutate(estimate = mu) %>%
mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>%
ggplot(aes(x= nitrate_concentration, y= estimate)) +
geom_point() +
geom_line(data=all_predsn_2, aes(x=nitrate_concentration.x, y=growth_rate, color = treatment), size = 1) +
facet_grid(treatment ~ ancestor_id) +
geom_hline(yintercept = 0.1, linetype = "dotted") +
ylab("Exponential growth rate (/day)") + xlab("Nitrate concentration (uM)")
monod_wide <- monod_fits %>%
select(population, term, estimate) %>%
spread(key = term, value = estimate)
m <- 0.1 ## set mortality rate, which we use in the rstar_solve
monod_curve_mortality <- function(nitrate_concentration, umax, ks){
res <- (umax* (nitrate_concentration / (ks+ nitrate_concentration))) - 0.1
res
}
Find R*
rstars <- monod_wide %>%
mutate(rstar = uniroot.all(function(x) monod_curve_mortality(x, umax, ks), c(0.0, 50))) %>% ## numerical
mutate(rstar_solve = ks*m/(umax-m)) ## analytical
rstars2 <- left_join(rstars, treatments, by = "population") %>%
distinct(population, ks, umax, .keep_all = TRUE)
rstars2 %>%
group_by(treatment) %>%
summarise_each(funs(mean, std.error), rstar_solve) %>%
ggplot(aes(x = reorder(treatment, mean), y = mean)) + geom_point() +
geom_errorbar(aes(ymin = mean - std.error, ymax = mean + std.error),width = 0.1) +
ylab("R* (umol N)") + xlab("Selection treatment") + geom_point(aes(x = reorder(treatment, rstar), y = rstar, color = ancestor_id), size = 2, data = rstars2, alpha = 0.5) +
scale_color_discrete(name = "Ancestor")